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Abstract. We use a transfer- matrix method to study the localization properties of vibrations 
in a 'mass and spring' model with simple cubic lattice structure. Disorder is applied as a box- 
distribution to the force-constants k of the springs. We obtain the reduced localization lengths 
Am from calculated Lyapunov exponents for different system widths to roughly locate the squared 
critical transition frequency ■ The data is finite-size scaled to acquire the squared critical transition 
frequency of u% = 12.54 ± 0.03 and a critical exponent of v — 1.55 ± 0.002. 



The disorder-induced metal-insulator transition (MIT) and the concept of Anderson localization 
[T] for electrons has been studied extensively for many years. During this time highly accurate results 
for the critical disorder and critical exponent have been obtained through numerous numerical 
methods, such as multifractal analysis [21 [3], energy level statistics [HE] and the transfer matrix 
method (TMM) [6j [7j. Although phonon localization has been a topic of research for a similar 
duration, where the first study of a disordered phonon is credited to Lifshitz et al. in 1954 [8], 
comparatively, less research has been conducted in the field. Analysis of the dynamical matrix 
through participation ratios [9j [10], level-spacing statistics [T2] and multifractal analysis [13] 
requires highly accurate numerical results on large system sizes which is not trivial and remains 
challenging even with todays computing resources. 

There are many features of interest in phononic research, such as the origins of the boson peak 
|lipi4j. the name given to the peak in the vibrational density of states (VDOS) g{oo) that is an excess 
contribution compared to the usual Debye behaviour g(u) oc uj 2 . The nature of the modes within 
the peak are of importance as they are strongly related to the mechanisms of thermal transport 
in the « 10 K temperature range [15]. Schumacher et al. [TTJ [15] successfully demonstrated that 
the boson peak consists of extended phonon modes. The locations of the boson peak, Ioffe-Regel 
transition between weak and strong elastic-scattering regimes and the localization-delocalization 
transition (LDT) are also of great interest in particular how these are related with each other. 
Scorpigno et al. [16J state that the link between the locations of the boson-peak and Ioffe-Regel 
transition is still debatable, whereas through level statistics Kantelhardt et al. [14] estimate the 
critical phonon frequency to be w c ~ 3cj pea k. Taraskin et al. [IT] declare that the links between 
vibrational localization, the boson peak, the Ioffe-Regel crossover and the zero-energy spectral 
singularity are yet to be fully established. 

Here, we take a simple cubic lattice structure made up of point masses connected by springs. 
We assume that the masses are displaced by only small amplitudes and apply the harmonic 



Figure 1. Amplitudes of vibration in a box of length N = 70 with periodic boundary conditions 
and evenly distributed force constant disorder with mean (k) = 1 and width Ak = 1. (a) Extended 
state lo 2 = 12, (b) critical state near the LDT transition to 2 = 12.5 and (c) localized state lo 2 = 13.03. 
Colours are an indication of the depth within the box that a particular amplitude is situated, the 
scale for which is included in the lower right hand side of the boxes. Sizes of internal boxes are 
determined from vibrational amplitudes at lattice site. Amplitudes less than the mean are not 
displayed. 

approximation. This gives rise to the classical model of vibration, where the displacements of 
the masses from their equilibrium positions are governed by the displacements of their nearest 
neighbours. Due to the symmetry of the lattice we can de-couple the spatial dimensions of the 
system into three identical problems and massively reduce the time of calculations. This is achieved 
by assuming that all springs have identical central and non-central force constants reducing the 
3x3 force constant matrices to scalars. Therefore, solving for just one spatial dimension achieves 
all three components of the vibrational amplitude for all lattice sites. 

We apply the above approximations and hence have to solve the standard equation of motion 
for a single lattice site with a plane wave solution to obtain 



where u is the frequency of vibration, m^ m ^ n and u\^ m ^ n are the mass and amplitude of vibration 
of the lattice site denoted by indices l,m,n. k is the force constant of the spring connecting two 
adjacent lattice sites, where the indices of the current lattice site I, m, n are omitted to simplify 
notation. The system of equations for all lattice sites is converted to matrix form by substitution 
into the generalised equation of motion 
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Figure 2. Reduced VDOS for 
evenly distributed spring disorder 
with mean (k) = 1 and width Ak = 
1. 150 bins are used and the data 
is normalized by total number of 
states. 



The above matrices [m] and [k] contain all the masses and force constants of the eigensystem 
and u is a vector of all vibrational amplitudes of the lattice sites. [m] _1 [k] is the dynamical matrix 
and is diagonalized to obtain the normal modes of vibration u with eigenvalues — w 2 . Although the 
vibrational amplitudes u are for only a single spatial dimension, they are also valid as the total 
magnitude of the amplitude. Recombining the individual amplitudes to obtain the total requires 
multiplication of single components by \^3, renormalization of the eigenvector removes these factors. 

Several methods can be used to analyse both the eigenvectors and eigenvalues obtained from 
diagonalization of the dynamical matrix to distinguish between extended and localized states. Just 
looking at the amplitude eigenvectors is a viable method of studying the localization lengths of 
the phonons, although it is not trivial due to the localization lengths diverging as the critical 
frequency is approached from the localized regime. Therefore in finite sized systems the localized 
states can appear extended as the localization length can be greater than the length of the system. 
We diagonalize the dynamical matrix for a box of size 70 3 with the inclusion of box distribution 
disorder with width AA; = 1 and mean (k) = 1 applied to the force constants of the springs. In figure 
QJa)-(c) we visualise three amplitude eigenvectors for squared frequencies of lo 2 = 12, w 2 = 12.5 and 
oj 2 = 13.03 respectively. Figure d{a) is that of an extended state as the amplitudes exhibit short 
range periodicity within the box and figure QJc) is that of a localized state where the amplitudes 
of vibration are confined to a small section of the box. Figure QJb) exhibits properties of both 
extended and localized states and is close to the transition. 

We calculate the reduced VDOS g(u)/u 2 for cubes with size 5 3 , 10 3 and 15 3 with 1360, 170 and 
50 random seeds respectively, resulting in roughly 170,000 states for each system size. The VDOS 
plot in figure [2] demonstrates the latter stages of the broadening of the Van Hove singularities due 
to disorder and the emergence of the characteristic low- frequency boson peak at Wp Ca k ~ 2. The 
frequencies co 2 = 12, 12.5 and 13.03 of the states plotted in figure CD are located close to the upper 
band edge of the VDOS, therefore we can safely assume that the LDT will be in this region. 

We introduce the TMM which calculates the localization length A of a quasi-one dimensional 
bar. For simplicity the bar geometry has an identical height and width M and in all cases the length 
of the bar L ^> M. L is not preset as the transfer matrix calculations continue along the length of 
L until the desired accuracy of the Lyapunov exponent is met, this convergence criterion is chosen 
before starting the calculation. 

To calculate the Lyapunov exponents, the system of equations ([I]) is converted to a form suitable 
for the TMM. This requires the equations to be re- arranged so that the amplitudes of a particular 
layer of lattice sites can be calculated solely from the amplitudes of the previous layer and the force 
constants of the surrounding springs. Computationally this is far less memory consuming than 
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Figure 3. Spring disorder AA: = 1.0: (a) reduced localization lengths for a range of frequencies 
and system sizes (12-20 in steps of 2). (b) Scaling function (solid line) and scaled data points for 
the expansion n r = 4, m r = 1, rii = 2, mj = 0. 



storing the whole system of masses and springs, as the TMM calculation only requires parameters 
for 2 layers of the system to be stored in memory at any one time. Equation ([1]) for a single lattice 
site is rearranged into a form where the amplitude of vibration of the nearest neighbour in layer 
I + 1 is calculated solely from parameters of sites in layers I and I — 1. 
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We introduce H\ = k^ min+ iu^ m , n+ i + fcz, m ,„-iu/,, m , n -i + ki,m+i,nn,m+l,n + k,m- 
as a collection of the in-plane contributions to the final amplitude and fc a n is the 
surrounding force constants for site l,m,n. We define £/>, Ui + i and ?7/_i as vectors 
the amplitudes of the constituent sites of layers I, I + 1 and / — 1 respectively. 
U\ = (^,1,1)^,1,2)^,2,1) • • • )^Z,M,Af)j where M is the width and height of the layer. 
([3]) is now expressed in the following matrix form, where Hi is a matrix containing 
contributions, and 1 are the zero and unit matrices, u 2 is single valued and all other 
are vectors. 
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The transfer matrix used to progress from one slice to the next is given as T ra and repeated 
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multiplication of this gives: Tjy = T n which according to Oseledec [IB] guarantees existence of 
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lim 



T N T N 



with eigenvalues e 7 . Where 7's are Lyapunov exponents and an 

estimation of localization length is given by the inverse of the minimum exponent X(M, u 2 ) = 1 /7 m in- 
The reduced (dimensionless) localization length may then be calculated as Aj\,/(cj 2 ) = A(M, u 2 )/M. 
We numerically calculate the reduced localization lengths for a range of squared frequencies and 
even system widths M between 12 and 20 to a convergence criterion of 0.1 percent of the variance. 
The reduced localization lengths Am as function of the squared frequency uj 2 are plotted in figure 
Ha). 



Table 1. Table of critical parameters obtained from finite size scaling for A A: = 1 on even numbered 
system widths between 12 and 20 and a range of squared frequencies u> 2 = 12.5 — 12.6. NDF is the 
number of degrees of freedom. We note that preliminary errors for v are too small as they have not 
been fully checked with respect to stability and robustness. 
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The reduced localization lengths from the TMM are finite size scaled by applying the procedure 
outlined in reference |19j . The finite size scaling (FSS) procedure obtains the correlation length £ 
for an infinite system from the reduced localization length Am {to 2 ) of finite sized systems by using 
the one-parameter scaling law Am = f(M/£) [20]. We assume that the LDT is characterised by 
a divergent correlation length, so that at a fixed disorder Ak, £(w 2 ) oc| uo 2 — lo 2 \~ v ', where v is 
the critical exponent and u 2 is the critical frequency squared. Highly accurate numerical studies 
of the Anderson model for electron localization have found the critical exponent v = 1.5 ± 0.1 
\2\\ [22| [T9] for box distribution disorder. Whereas in the phononic case Monthus and Garel [23] 
have demonstrated that their participation ratio data for high disorder at an LDT collapsed fairly 
well onto their scaling function using a critical exponent u = 1.57 taken from reference [21]. Akita 
and Ohtsuki [25J find a critical exponent of v ~ 1.2 ± 0.2 by finite size scaling TMM data at an 
LDT for box distribution spring disorder Ak = 1.8 with a convergence criterion of 2 percent. Table 
Q] contains the critical parameters obtained from FSS for particular orders of Taylor expansions of 
the relevant and irrelevant scaling variables n t , rii, m r and m\ listed in columns 1-4. Figure E^b) 
is the scaling function obtained from the highest order stable fit given in tabled] (n T = 4, m = 1, 
m r = 2, rri[=0). The mean critical exponent (v) = 1.55 ± 0.002 which is remarkably close to the 
MIT critical exponent. 

From FSS we achieve a weighted mean critical squared transition frequency of lo 2 = 12.54 ±0.03 
which as we can see from figure [2] occurs in the tail of the VDOS and close to the example state of 
squared frequency lo 2 = 12.5 in Figure QJb) that exhibits properties of both extended and localized 
states. It is generally assumed that the LDT shares the same exponent as the MIT, where the 
current accepted value is f « 1.5±0.1. Although Akita and Ohtsuki [25 j have previously obtained 
a much lower value v ~ 1.2 ± 0.2, thought to be a consequence of the transition occurring close 
to the upper band edge [26], we have achieved a striking resemblance between critical exponents 
between the MIT and LDT. Therefore we assume that the phononic model will be of the same 
orthogonal universality class as the Anderson model. 
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